library(Matching)
library(foreign)
library(rbounds)
rm(list=ls())

setwd("/Users/guy/Dropbox/Dissertation/Voting Rule/QJPS/replication files/")
source("HodgesLehman.R")

# Define Outcome variables
y_names <- c(
  friendA = "Network centrality sumamry index",               
  monitoringA = "Monitoring summary index",             
  cooperationA = "Coopeartion summary index", 
  responsivenessA = "Responsiveness summary index",  
  z_imsf33_memb = "Manager is accountable (M)",
  z_imsf32_d2_memb  = "Manager is monitored (M)",  
  z_sf222_reps = "Manager is monitored (R)",
  z_isd15_memb = "General assembly: members' attendance (M)",
  z_isd15_reps = "General assembly: reps' attendance (R)",  
  z_irsf21_reps = "Internal Auditing (R)",
  z_irsf20b_reps = "External auditing (R)  ",
  z_irsf10_reps = "Representatives influence over decisions (R)",
  z_imse09b_memb = "Members dries coffee on tarps",
  z_isb13_memb ="Members planted sidling in past 12 months",   
  z_imbulk_a_memb= "Members' share of coffee bulked",  
  z_contribution_memb = "Members' contribution: commitment experiment",  
  z_isd20_memb = "Members paid annual dues",  
  z_isd18_memb = "Members paid joining fees",
  z_isd28_memb = "Members agree to increase commission",
  z_imsf28_d4_memb = "Manager is very transparent (M)",
  z_imse11_memb = "Members warned: bad agricultural practices (M)",  
  z_ss_warned2_memb  = "Members warned: side selling (M)",
  z_irsf19b_reps =  "Receipts are given to members (R)",
  z_imsf30b_memb = "Members obtain receipts (M)",
  isd12b_memb = "pcnt meetings attended past season (MS)",                
  isd06b_memb = "pcnt trainings attended past season (MS)"              
 )

DC <- read.dta("GenMatch.dta")
attach(DC)
head(DC)

# Create matched pairs
set.seed(20111210)
Y  <- DC$monitoringA
Tr <- DC$z
X <- cbind(DC$agedc, DC$npos_before, DC$dsbq8b)
BalanceMatrix <- cbind(DC$agedc + I(DC$agedc^2) + DC$npos_before +  I(DC$npos_before^2) + DC$dsbq8b +  I(DC$dsbq8b^2) + I(DC$agedc*DC$npos_before) + I(DC$agedc*DC$dsbq8b))
gen1 <- GenMatch(Tr = Tr, X = X, BalanceMatrix = BalanceMatrix, pop.size = 1000, replace=FALSE, paired=TRUE)

# obtain balance statistic
mgen1 <- Match(Y = Y, Tr = Tr, X = X, Weight.matrix = gen1, replace=FALSE)
summary(mgen1)
mgen1$index.treat
psens.test <- psens(mgen1, Gamma=2, GammaInc=.1) # conduct sensitvity analysis of p-values
print(psens.test)
hsens.test<- hlsens(mgen1, Gamma=2, GammaInc=.1, .1) # Rosenbaum Sensitivity Test for Hodges-Lehmann Point Estimat
print(hsens.test)
MatchBalance(Tr ~ DC$agedc + I(DC$agedc^2) + DC$npos_before +  I(DC$npos_before^2) + DC$dsbq8b +  I(DC$dsbq8b^2) + I(DC$agedc*DC$npos_before) + I(DC$agedc*DC$dsbq8b), match.out =  mgen1, nboots = 1000, data=DC)

DC_ctrl <- DC[mgen1$index.control,]
DC_trt  <- DC[mgen1$index.treated,]

# Mann-Whitney test for IV 
set.seed(313457)
foo <- function(y_name, conf.level,...) {
  DC <- DC[!is.na(DC[[y_name]]),] 
  wtiv <- wilcox.test(x = DC_trt[[y_name]], y = DC_ctrl[[y_name]], paired = TRUE, conf.int=FALSE, exact=FALSE)

# define IV function
  complete <- na.omit(cbind(DC_trt[[y_name]], 
                            DC_ctrl[[y_name]], 
                            DC_trt$vrule, 
                            DC_ctrl$vrule, 
                            DC_trt$z, 
                            DC_ctrl$z)) 
  
  HL <-  HodgesLehman(Rt = complete[,1],
                      Rc = complete[,2],
                      Dt = complete[,3],
                      Dc = complete[,4],
                      Zt = complete[,5],
                      Zc = complete[,6])
  
  out <- 
    iv_sens(Rt = complete[,1],
            Rc = complete[,2],
            Dt = complete[,3],
            Dc = complete[,4],  Gamma = 4, GammaInc = .1) 
  
  out$wt_p.value <- wtiv[["p.value"]]
  out$matched_pairs <- nrow(DC_trt)
  
  DC_ctrl$id <- 1:nrow(DC_ctrl)
  DC_trt$id  <- 1:nrow(DC_trt)
  
  out$control <- mgen1$index.control
  out$treated <- mgen1$index.treated
  out$pairs <- rbind(DC_ctrl, DC_trt)
  out$pairs <- out$pairs[order(out$pairs$id),]
  
  out$HL <- HL$HL
  out$pval <- HL$pval
  out$Interval <- HL$Interval 
  return(out)
}

# define ITT function
bar <- function(y_name, conf.level,...) {
  diff<-as.vector(DC_trt[[y_name]] - DC_ctrl[[y_name]])
  diff <- diff[!is.na(diff)]
  if(length(unique(diff)) < 2) return(NA)
  wt <- try(wilcox.test(diff, conf.int=TRUE), silent = FALSE)
  if(class(wt)[1] == "htest") return(wt)
  else return(NA)
}

# execute IV and ITT functions
dim(DC)
include <- 6:32
new <- colnames(DC[,include])
iv_025 <- sapply(new, FUN = foo, alpha = .025, conf.level=0.95, Beta = 2, simplify = FALSE)
iv_025

new <- colnames(DC[,include])
itt_025 <- sapply(new, FUN = bar, alpha = .025, conf.level=0.95, Beta = 2, simplify = FALSE)
itt_025 

# create matrix of IV and ITT results
outcomes <- names(iv_025)
mat1 <- matrix(NA,length(outcomes),8)
rownames(mat1) <- outcomes
for (i in 1:length(outcomes)) {
  mat1[i,1] <- iv_025[[i]]$pval
  mat1[i,2] <- iv_025[[i]]$HL
  mat1[i,3:4] <- c(iv_025[[i]]$Interval[1,1], iv_025[[i]]$Interval[1,2])
  if(length(itt_025[[i]]) == 1) next
  mat1[i,5] <- itt_025[[i]]$p.value
  mat1[i,6] <- itt_025[[i]]$estimate
  mat1[i,7:8] <- itt_025[[i]]$conf.int
}

mat1[is.infinite(mat1)] <- NA_real_
colnames(mat1) <- c(paste("IV:", c("pval", "estimate", "lower", "upper")),
                    paste("ITT:", c("pval", "estimate", "lower", "upper")))

# export results to excel for graphing
write.csv(mat1,"GenMatchingEstimation.csv")
write.csv(cbind(trt = mgen1$index.treated, ctrl = mgen1$index.control), "pairs.csv")
cat("success")

# Sensitivity analysis
sensitivityA <- iv_025[["cooperationA"]]
sensitivityB <- iv_025[["responsivenessA"]]
sensitivityC <- iv_025[["monitoringA"]]
sensitivityD <- iv_025[["z_group_bias"]]

par(mfrow=c(2,2))
plot(sensitivityA$bounds[,-2],
     xlab = "Gamma", ylab = "p-value", 
     pch = 20, font.lab = 2, ylim = c(0,0.25), main="Cooperation Index", 
     abline(h = 0.05, lty = 2))

plot(sensitivityB$bounds[,-2],
     xlab = "Gamma", ylab = "p-value", 
     pch = 20, font.lab = 2, ylim = c(0,0.25), main="Responsiveness Index", 
     abline(h = 0.05, lty = 2))

plot(sensitivityC$bounds[,-2],
     xlab = "Gamma", ylab = "p-value", 
     pch = 20, font.lab = 2, ylim = c(0,0.25), main="Monitoring Index", 
     abline(h = 0.05, lty = 2))

plot(sensitivityD$bounds[,-2],
     xlab = "Gamma", ylab = "p-value", 
     pch = 20, font.lab = 2, ylim = c(0,0.25), main="Member-Stranger (DG)", 
     abline(h = 0.05, lty = 2))
dev.off()
